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We consider a mesoscopic region coupled to two leads under the influence of 
external time-dependent voltages. The time dependence is coupled to source 
and drain contacts, the gates controlling the tunnel-barrier heights, or to 
the gates which define the mesoscopic region. We derive, with the Keldysh 
nonequilibrium Green function technique, a formal expression for the fully 
nonlinear, time-dependent current through the system. The analysis admits 
arbitrary interactions in the mesoscopic region, but the leads are treated as 
noninteracting. For proportionate coupling to the leads, the time-averaged 
current is simply the integral between the chemical potentials of the time- 
averaged density of states, weighted by the coupling to the leads, in close 
analogy to the time-independent result of Meir and Wingreen (Phys. Rev. 



Lett. 68, 2512 (1992)). Analytical and numerical results for the exactly 
solvable non- interacting resonant-tunneling system are presented. Due to the 
coherence between the leads and the resonant site, the current does not follow 
the driving signal adiabatically: a "ringing" current is found as a response to a 
voltage pulse, and a complex time-dependence results in the case of harmonic 
driving voltages. We also establish a connection to recent linear-response 
calculations, and to earlier studies of electron-phonon scattering effects in 
resonant tunneling. 

73.20.Dx, 73.40.Ei, 73.40.Gk, 73.50.Fq 



2 



Typeset using REVT^X 



I. INTRODUCTION 



The hallmark of mesoscopic phenomena is the phase coherence of the charge carriers, 
which is maintained over a significant part of the transport process. The interference effects 
resulting from this phase coherence are reflected in a number of experimentally measur- 
able properties. For example, phase coherence is central to the Aharonov-Bohm effect, [I[ 
Universal conductance fluctuations, |T] and weak localization, and can be affected by 
external controls such as temperature or magnetic field. The study of stationary mesoscopic 
physics is now a mature field, and in this work we focus on an alternative way of affecting 
the phase coherence: external time- dependent perturbations. The interplay of external time 
dependence and phase coherence can be phenomenologically understood as follows. If the 
single-particle energies acquire a time dependence, then the wave functions have an extra 
phase factor, ip ~ exp(— i J t dt'e(t')). For a uniform system such an overall phase factor is 
of no consequence. However, if the external time dependence is different in different parts 
of the system, and the particles can move between these regions (without being 'dephased' 
by inelastic collisions), the phase difference becomes important. 

The interest in time-dependent mesoscopic phenomena stems from recent progress in 
several experimental techniques. || Time dependence is a central ingredient in many different 
experiments, of which we mention the following: (i) Single- electron pumps and turnstiles. 
Here time-modified gate signals move electrons one by one through a quantum dot, leading to 
a current which is proportional to the frequency of the external signal. These structures have 
considerable importance as current standards. The Coulombic repulsion of the carriers in the 
central region is crucial to the operational principle of these devices, and underlines the fact 
that extra care must be paid to interactions when considering time-dependent transport in 
mesoscopic systems, (ii) ac response and transients in resonant-tunneling devices. Resonant 
tunneling devices (RTD) have a number of applications as high-frequency amplifiers or 
detectors. For the device engineer a natural approach would be to model these circuit 
elements with resistors, capacitances, and inductors. The question then arises as to what, if 
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any, are the appropriate 'quantum' capacitances and inductances one should ascribe to these 
devices. Answering this question requires the use of time-dependent quantum-transport 
theory, (iii) Interaction with laser fields. Ultrashort laser pulses allow the study of short- 
time dynamics of charge carriers. Here again, coherence and time dependence combine with 
the necessity of treating interactions. 

A rigorous discussion of transport in an interacting mesoscopic system requires a formal- 
ism which is capable of including explicitly the interactions. Obvious candidates for such a 
theoretical tool are various techniques based on Green functions. Since many problems of 
interest involve systems far from equilibrium, we cannot use linear-response methods, such 
as those based on the Kubo formula, but must use an approach capable of addressing the 
full nonequilibrium situation. The nonequilibrum Green function techniques, as developed 
about thirty years ago by Kadanoff and Baym, H and by Keldysh, have during the re- 
cent years gained increasing attention in the analysis of transport phenomena in mesoscopic 
semiconductors systems. 0j In particular, the steady state situation, has been addressed 
by a large number of papers [§HTB| Among the central results obtained in these papers is 
that that under certain conditions (to be discussed below) a Landauer type conductance 
formula WM can be derived. This is quite appealing in view of the wide spread success of 



conductance formulas in the analysis of transport in mesoscopic systems. 

Considerably fewer studies have been reported where an explicit time dependence is an 
essential feature. We are aware of an early paper in surface physics, |L5| but only in the 



recent past have groups working in mesoscopic physics addressed this problem. [|16| -|2T|j The 
work reported in this paper continues along these lines: we give the full details and expand 



on our short communication. 18 



Our main formal result from the nonequilibrium Green function approach is a general 
expression for the time-dependent current flowing from non-interacting leads to an inter- 
acting region. As we will discuss in Section II, the time dependence enters through the 
self-consistent parameters defining the model. We show that under certain restrictions, to 
be specified below, a Landauer-like formula can be obtained for the time-averaged current. 



To illustrate the utility of our approach we give results for an exactly solvable non-interacting 
case, which displays an interesting, and experimentally measurable, nonadiabatic behavior. 
We also establish a link between the present formulation and recently published results for 
linear-response and electron-phonon interactions, obtained by other techniques. 

The paper is organized as follows. We examine in Section II the range of experimental 
parameters in which we expect our theoretical formulation to be valid. In Section III we 
briefly review the physics behind the nonequilibrium Green function technique of Keldysh, 
and Baym and Kadanoff, which is our main theoretical tool, and then introduce the specific 
model Hamiltonians used in this work. We derive the central formal results for the time- 
dependent current in Section IV. We also derive, under special restrictions, a Landauer- 
like formula for the average current. In Section V we apply the general formulae to an 
explicitly solvable resonant-tunneling model. Both analytical and numerical results are 
presented. We also show that the linear ac-response results of Fu and Dudley p2 
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contained as a special case of the exact results of this section. In Section VI we illustrate the 
utility of our formulation by presenting a much simplified derivation of Wingreen et al's 
results on resonant tunneling in the presence of electron-phonon interactions. Appendix A 
summarizes some of the central technical properties of the Keldysh technique: we state the 
definitions, give the basic equations, and provide the analytic continuation rules employed 
below. In Appendices B and C we present proofs for certain statements made in the main 
text, and, finally, in Appendix D we describe some transformations which facilitate numerical 
evaluation of the time- dependent current. 



II. APPLICABILITY TO EXPERIMENTS 

A central question one must address is: under which conditions are the non-equilibrium 
techniques, applied successfully to the steady-state problem, transferrable to time-dependent 
situations, such as the experiments mentioned above? 

The time-dependent problem has to be formulated carefully, particularly with respect 
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to the leads. It is essential to a Landauer type of approach, that the electrons in the leads 
be non-interacting. In practice, however, the electrons in the leads near the mesoscopic 
region contribute to the self-consistent potential. We approach this problem by dividing the 
transport physics in two steps: (i) the self-consistent determination of charge pile-up 
and depletion in the contacts, the resulting barrier heights, and single-particle energies in 
the interacting region, and (ii) transport in a system defined by these self-consistent pa- 
rameters. Step (i) requires a capacitance calculation for each specific geometry, [|4[] and we 
do not address it in this paper. Instead, we assume the results of (i) as time- dependent 
input parameters and give a full treatment of the transport through the mesoscopic region 
(ii). In practice, the interactions in the leads are absorbed into a time-dependent potential 
and from then on the electrons in the leads are treated as non-interacting. This means that 
when relating our results to actual experiments some care must be exercised. Specifically, we 
calculate only the current flowing into the mesoscopic region, while the total time-dependent 
current measured in the contacts includes contributions from charge flowing in and out of 
accumulation and depletion regions in the leads. In the time- averaged (dc) current, how- 
ever, these capacitive contributions vanish and the corresponding time-averaged theoretical 
formulae, such as Eq. (]27|) , are directly relevant to experiment. It should be noted, though, 
that these capacitive currents may influence the effective time-dependent parameters in step 
(i) above. 

Let us next estimate the frequency limits that restrict the validity of our approach. Two 
criteria must be satisfied. First, the driving frequency must be sufficiently slow that the 
applied bias is dropped entirely across the tunneling structure. When a bias is applied to a 
sample, the electric field in the leads can only be screened if the driving frequency is smaller 
than the plasma frequency, which is tens of THz in typical doped semiconductor samples. 
For signals slower than this, the bias is established entirely across the tunneling structure by 
accumulation and depletion of charge near the barriers. The unscreened Coulomb interaction 
between net excess charge is quite strong, and hence the bias across a tunneling structure 
is caused by a relatively small excess of charge in accumulation and depletion layers. The 
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formation of these layers then causes a rigid shift (see Eq.(§) below) of the bottom of the 
conduction band deeper in the leads, which is the origin of the rigid shift of energy levels in 
our treatment of a time-dependent bias. 

The second frequency limit on our approach is that the build-up of electrons required 
for the formation of the accumulation and depletion layers must not significantly disrupt 
the coherent transport of electrons incident from the leads. One way to quantify this is to 
ask - what is the probability that an electron incident from the leads participates in the 
build-up of charge associated with a time- dependent bias? This probability will be the ratio 
of the net current density flowing into the accumulation region to the total incident flux of 
electrons. For a three dimensional double-barrier resonant-tunneling structure (see Fig. |]) 
the ac-current charging the accumulation layer is J™ s = 2irvCV rras / A, where v is the driving 
frequency, C is the capacitance, V rms is the applied bias, and A is the area. In comparison, 
the total incident flux is L mc = 3/8envp. Using the parameters appropriate for a typical 
experiment (we use that of Brown et al. we find that up to 10 THz the probability of 
an electron participating in the charge build-up is only 1%. Summarizing, these estimates 
indicate that our approach should be accurate up to frequencies of tens of THz, which are 
large by present experimental standards [Ref. ?], and consequently the analysis presented 
in what follows should be valid for most experimental situations. 

III. THEORETICAL TOOLS, AND THE MODEL 

A. Baym-Kadanoff-Keldysh nonequilibrium techniques 

Here we wish to outline the physical background behind the Keldysh formulation, and 
in particular its connection to tunneling physics. Readers interested in technical details 
should consult any of the many available review articles, such as Refs. [ [p5|-|27|1. The basic 
difference between construction of equilibrium and nonequilibrium perturbation schemes is 
that in nonequilibrium one cannot assume that the system returns to its ground state (or 
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a thermodynamic equilibrium state at finite temperatures) as t — > +00. Irreversible effects 
break the symmetry between t = —00 and t = +00, and this symmetry is heavily exploited 
in the derivation of the equilibrium perturbation expansion. In nonequilibrium situations 
one can circumvent this problem by allowing the system to evolve from — 00 to the moment of 
interest (for definiteness, let us call this instant to), and then continues the time evolvement 
from t = to back to t — —00. |2£| (When dealing with quantities which depend on two time 
variables, such as Green functions, the time evolution must be continued to the later time.) 
The advantage of this procedure is that all expectation values are defined with respect to a 
well defined state, i.e. the state in which the system was prepared in the remote past. The 
price is that one must treat the two time branches on an equal footing (See Fig. 0). 

A typical object of interest would be a two time Green function (see Appendix A); the 
two times can be located on either of the two branches of the complex time path (e.g., 
r and r' in Fig. |2|). One is thus led to consider 2x2 Green function matrices, and the 
various terms in the perturbation theory can be evaluated by matrix multiplication. Since 
the internal time-integrations run over the complex time path, a method of book-keeping 
for the time-labels is required, and there are various ways of doing this. In the present work 
we employ a version of the Keldysh technique. 

In the context of tunneling problems the time-independent Keldysh formalism works as 
follows. In the remote past the contacts (i.e. the left and right lead) and the central region are 
decoupled, and each region is in thermal equilibrium. The equilibrium distribution functions 
for the three regions are characterized by their respective chemical potentials; these do not 
have to coincide nor are the differences between the chemical potentials necessarily small. 
The couplings between the different regions are then established and treated as perturbations 
via the standard techniques of perturbation theory, albeit on the two-branch time contour. 
It is important to notice that the couplings do not have to be small, e.g. with respect level 
spacings or k^T, and typically must be treated to all orders. 

The time-dependent case can be treated similarly. Before the the couplings between 
the various regions are turned on, the single-particle energies acquire rigid time- dependent 

8 



shifts, which, in the case of the non-interacting contacts, translate into extra phase factors 
for the propagators (but not in changes in occupations). The perturbation theory with 
respect to the couplings has the same diagrammatic structure as in the stationary case. The 
calculations, of course, become more complicated because of the broken time translational 
invariance. 

B. Model Hamiltonian 

We split the total Hamiltonian in three pieces: H = H c + H T + H ccn , where H c describes 
the contacts, is the tunneling coupling between contacts and the interacting region, and 
H cen models the interacting central region, respectively. Below we discuss each of these 
terms. 

1. Contacts, H c 

Guided by the typical experimental geometry in which the leads rapidly broaden into 
metallic contacts, we view electrons in the leads as non-interacting except for an overall 
self-consistent potential. Physically, applying a time-dependent bias (electrostatic-potential 
difference) between the source and drain contacts means that the single-particle energies 
become time-dependent: e° Q — > €k a (t) = + A a (t) (here a labels the channel in the 
left (L) or right (R) lead). The occupation of each state ka, however, remains unchanged. 
The occupation, for each contact, is determined by an equilibrium distribution function 
established in the distant past, before the time dependence or tunneling matrix elements 
are turned on. Thus, the contact Hamiltonian is 

H c = e ka(t)cl a c ka , (1) 

k,adL,R 

and the exact time-dependent Green functions in the leads for the uncoupled system are 
g< a (t,t')=i(cl a (t')c ka (t)) 
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rt 

= if (4a) exp [ - i / dhekaiti)) 
Jt' 

gZ{t,t>)=TiO{±tTt'){{cUt)A«{m 

= T id(±t t exp [ - i f f dhekvih)] . (2) 

One should note that our model for g K differs from the choice made in the recent study of 



Chen and Ting: [If]] these authors allow the electrochemical potential in the distribution 
function / to vary with time: fii — = e[V + U(t)], where U(t) is the ac signal. This 
assumption implies that the total number of electrons in the contacts varies with time. 
This behavior is inconsistent of what happens in real devices: it is only the relatively small 
number of electrons in the accumulation/ depletion layers that is time- dependent. In addition 
to the unphysical charge pile-up in the contacts, the model of Chen and Ting leads to an 
instantaneous loss of phase-coherence in the contacts, and hence does not display any of the 
interesting interference phenomena predicted by our phase-conserving model. 



2. Coupling between leads and central region, Hf 

The coupling between the leads and the central (interacting) region can be modified 
with time-dependent gate voltages, as is the case in single-electron pumps. The precise 
functional form of the time- dependence is determined by the detailed geometry and by the 
self-consistent response of charge in the contacts to external driving. We assume that these 
parameters are known, and simply write 

H T = £ [V ka , n (t)cl a d n + h.c.}. (3) 

k,a£L,R 
n 

Here {d^} and {d„} form a complete orthonormal set of single-electron creation and anni- 
hilation operators in the interacting region. 
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3. The central region Hamiltonian H cen 

The form chosen for H cen in the central interacting region depends on geometry and on 
the physical behavior being investigated. Our results relating the current to local properties, 
such as densities of states and Green functions, are valid generally. To make the results more 
concrete, we will discuss two particular examples in detail. In the first, the central region is 
taken to consist of noninteracting, but time-dependent levels, 

H cen = J2^m(t)dl n d m (4) 

m 

Here d] n (d m ) creates (destroys) an electron in state m. The choice (|]) represents a simple 
model for time-dependent resonant tunneling. Below we shall present general results for an 
arbitrary number of levels, and analyze the case of a single level in detail. The latter is 
interesting both as an exactly solvable example, and for predictions of coherence effects in 
time-dependent experiments. 

The second example we will discuss is resonant tunneling with electron-phonon interac- 
tion, 

Hf~^ = 6 dtd + dtd £ M q [a q + a_ q ] , (5) 

q 

In the above, the first term represents a single site, while the second term represents in- 
teraction of an electron on the site with phonons: a q (a q ) creates (destroys) a phonon in 
mode q, and M q is the interaction matrix element. The full Hamiltonian of the system 
must also include the free-phonon contribution = J2q ^ q a q a q . This example, while 
not exactly solvable, is helpful to show how interactions influence the current. Furthermore, 



we can directly compare to previous time-independent results []23[ using (y) to demonstrate 
the power of the present formalism. 

IV. TIME-DEPENDENT CURRENT AND KELDYSH GREEN FUNCTIONS 
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A. General expression for the current 



The current from the left contact through left barrier to the central region can be calcu- 
lated from the time evolution of the occupation number operator of the left contact: 

ip 

J L (t) = -e(N L ) = --([H,N L ]), (6) 

where Nl = J2k,aeL c la c ka and H = H c + H T + H cen . Since H c and H cen commute with N L , 
one readily finds 

IP x 

J l = -r E [VkaMM - ^„<dic fa >] . (7) 

^ k,a£L 

71 

Now define two Green functions: 

G< ka (t,t')=z(cl a (t')d n (t)) (8) 

G< ajn (t,t')=i(di(t')c ka (t))- (9) 
r i * 

Using G^ an (t,t) = — G^ ka (t,t) , and inserting the time labels, the current can be ex- 
pressed as 

j L (t) = ^ Re { £ V ka>n (t)G< ka (t,t)} . (10) 

71 

One next needs an expression for G< ka (t,t'). For the present case, with non-interacting 
leads, a general relation for the contour-ordered Green function G n ^ a ( T , T ') can be derived 
rather easily (either with the equation-of-motion technique, or by a direct expansion of the 
S- matrix; the details are given in Appendix B), and the result is 

G nM {T,T') = E / dT l G nrn{T'^ 1 )V k \ m {T l )g ka {T l ,T') . (11) 

m 

Here G nm (r, n) is the contour-ordered Green function for the central region, and the r- 
variables are now defined on the contour of Fig. 0. Note that the time-dependence of the 
tunneling matrix elements and single-particle energies has broken the time-translational 
invariance. The analytic continuation rules ( |A3| ) of Appendix A can now be applied, and 
we find 
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m J 

+G< m (t,t 1 )g a ka (t l ,t')] , (12) 



where the Green functions g <,a for the leads are defined in (|2]) above. Combining (P) , (|i~0D 
and flT2|), yields 

2e 



= ~Im{ £ V ka , n (t) I dtjk***"™ 

11 k, a &L J -°° 

*V: ajm {h)[Gl m {tM)h{e ka ) + G< m (t,t 1 )}} . (13) 

The discrete sum over k in J2ka can be expressed in terms of densities of states in the leads: 
/ dep a (e). Then it is useful to define 



[r L Mi,t)] mn = 2tt £ p Q (6)y a , n (6,t)^%(e,ti) 

a£L 

x exp[i f dt 2 A a {e,t 2 )} , (14) 



where V ka>n = V^ )TV (ejfe). In terms of this generalized line-width function ([14]), the general 
expression for the current is 

n J-oo J 2n L 
x[G<(* ) * 1 ) + / i (c)G'-(*, «!)]}. (15) 

Here the bold-face notation indicates that the level-width function T and the central-region 
Green functions G <,r are matrices in the central-region indeces m, n. An analogous formula 
applies for Jn(t), the current flowing into the central region through the right barrier. 

This is the central formal result of this work, and the remainder of this paper is devoted to 
the analysis and evaluation of Eq . (|i~5|) . The current is expressed in terms of local quantities: 
Green functions of the central region. The first term in Eq.(|l5|), which is proportional 
to the lesser function G < , suggests an interpretation as the out-tunneling rate (recalling 
lmG < (t,t) = N(t)). Likewise, the second term, which is proportional to the occupation 
in the leads and to the density of states in the central region, can be associated to the 
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in-tunneling rate. However, one should bear in mind that all Green functions in Eq. (|15D 
are to be calculated in the presence of tunneling. Thus, G K will depend on the occupation 
in the leads. Furthermore, in the presence of interactions G r will depend on the central 
region occupation. Consequently, the current can be a non-linear function of the occupation 
factors. This issue has recently been discussed also by other authors. 



B. Time-independent case 

1. General expression 

In the time- independent limit the line- width function simplifies: T(e,ti,t) — > T(e), and 
the ti-integrals in Eq.fll^) can be performed: 

C dtJ ^lmTi\e- l ^- t) r L (e)G < (t-t 1 )\ 

= 4/^ Tr{rL(e)G<(e)} ' (16) 

and 

/I dtl I ^ l ^{^' t)vL ^)Ue)G r {t - h)} 

= -\J ^{r L (e)h(e)[G r (e) - G a (e)}} . (17) 

When these expressions are substituted to Eq.(|T5D, the current from left (right) contact to 
central region becomes 

^>4/^{ riM M< G< < £ > 

+/ iW ( e )[G'( £ )-G«( e )])}. (18) 

In steady state, the current will be uniform, so that J = Jl = —Jr , and one can symmetrize 
the current: J = (Jl + Jl)/2 = (Jl — Jr)/^- Using Eq . (|18|) leads to the general expression 
for the d.c. current: 

+ [h(e)T L (e) - f R (e)T R (e)}[G r (e) - G a (e)}} . (19) 

14 



This result was reported in Ref. @, and applied to the out-of-equilibrium Anderson impurity 
problem. 

2. Proportionate coupling 

If the left and right line-width functions are proportional to each other, i.e. T L (e) = 
XT R (e), further simplification can be achieved. We observe that the current can be written 
as J = xJl — (1 — x)Jr, which gives, using Eq.([H|), 

J=jj^Tr{T R (e)[(Xx-(l-x))G<(e) 

+ (Xxf L - (1 - x)f R )(G r (e) - G a )(e)]} . (20) 

The arbitrary parameter x is now fixed so that the first term vanishes, i.e x — 1/(1 + A), 
which results in 

xT '1 r^)Tr«( e ) (Gr(6> " G ° (e)> ^ (21) 

The ratio is well-defined because the T-matrices are proportional. The difference between 
the retarded and advanced Green functions is essentially the density of states. Despite of 
the apparent similarity of (^T|) to the Landauer formula, it is important to bear in mind that 
in general there is no immediate connection between the quantity Tr|(r L (e)r i? (e)/[r i (e) + 
r i? (e)])[G r (e) — G a (e)]|, and the transmission coefficient. In particular, when inelastic 
scattering is present, we do not believe that such a connection exists. In Section V, where 
we analyze a non-interacting central region, a connection with the transmission coefficient 
can be established. Further, in the next section we shall see how an analogous result can be 
derived for the average of the time-dependent current. 



15 



C. Average current 



In analogy with the previous subsection, where we found a compact expression for the 
current for the case of proportionate coupling, the time-dependent case allows further sim- 
plification, if assumptions are made on the line-width functions. In this case, we assume a 
generalized proportionality condition: 

r L (e,t 1 ,t) = XT R (e,t 1 ,t). (22) 

One should note that in general this condition can be satisfied only if A%(t) = A R (t) = 
A(t). However, in the Wide-Band Limit (WBL), to be considered in detail below, the 
time-variations of the energies in the leads do not have to be equal. 

We next consider the occupation of the central region N(t) = Sm(dJ n (t)d m (t)) and apply 
the continuity equation: 

dN(t) 
6 dt 

which allows one to write for arbitrary x: 



Mt) + J L (t) , (23) 



J L (t) = xJ L (t) + (1 - " Mt)] • (24) 



Choosing now x = 1/(1 + A) leads to 



Mt) ' A 



e— -rlmTrj / dt\ [ — 

dt n l J-oo J 2tt 



x 



1 + A 

e-Hti-t^R^^^r^^f^ _ fR{e)] j 



(25) 



The time-average of a time-dependent object F(t) is defined by 

= £l dtF ® • (26) 



If F(t) is a periodic function of time, it is sufficient to average over the period. Upon time- 
averaging, the first term in Eq. fl25D vanishes, (dN/dt) — > 0, because the occupation N(t) is 
finite for all T. The expression for the time-averaged current further simplifies if one can 
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factorize the energy- and time-dependence of the tunneling coupling, Vfc Q)TV (t) = u (t) V Qjn (efc). 
We then obtain 

XIm ^{ r^+ r r"t) < " ( ' )A(e ' f)> ^ (27) 

where 

A(e,t) = / dtiM(ti)G r (t,t!) 

x exp[ze(t - t x ) + i /* dt 2 A(t 2 )] . (28) 

Due to Eq.(|22|) we do not have to distinguish between L/R in the definition of A(e, £); below 
we shall encounter situations where this distinction is necessary. 

The expression (|27|) is of the Landauer type: it expresses the current as an integral over 
a weighted density of states times the difference of the two contact occupation factors. It is 
valid for arbitrary interactions in the central region, but it was derived with the somewhat 
restrictive assumption of proportional couplings to the leads. 



V. NON-INTERACTING RESONANT-LEVEL MODEL 

A. General formulation 

In the no n- interacting case the Hamiltonian is H = H c + + H ccn , where H cen = 
Y^n e ndnd n . Following standard analysis (an analogous calculation is also carried out in 
Appendix B), one can derive the Dyson equation for the retarded Green function, 

G r (t,t') =g r (t,t') 

+ jdhj d^ 2 g'•(^,^ 1 )s ^ ■(^ 1 ,^ 2 )G ^ •(^ 2 ,^ , ) , (29) 



where 



£; n ,(t l5 t 2 )= E yL^UhMVkaA^) , (30) 

ka£L,R 
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and g r ka is given by Eq.(g). From ( |A4j) the Keldysh equation for G < is p2 



G<(t,t') = J dhj rft 2 G r (t,t 1 )S < (t 1 ,t 2 )G a (t 2 ,t' 



^i e -ie( tl -t 2 ) 



f L/R {e)T L l R {eM,h 



G a (t 2 ,t') . 



(31) 



Provided that the Dyson equation for the retarded Green function can be solved, Eq.(|3T|) 
together with the current expression Eq. flTol) provides the complete solution to the nonin- 
teracting resonant-level model. Below we examine special cases where analytic progress can 
be made. 



B. Time-independent case 

In the time-independent case the time-translational invariance is restored, and it is ad- 
vantageous to go over to energy variables: 

G<(e) = G r (e)12<(e)G a (e) . (32) 

In general the Dyson equation for the retarded Green function requires matrix inversion. In 
the case of a single level, the scalar equations can be readily solved. The retarded (advanced) 
self-energy is 

S7-"(e)= £ lVkal [. = A(e) =f k(e) , (33) 
telL,B e - <*a ± wj 2 

where the real and imaginary parts contain 'left' and 'right' contributions: A(e) = A L (e) + 
A R (e) and T(e) = r L (e) + T R (e). The lesser self-energy is 

£<( e )= E \Vka\ 2 9^(e) 

kaeL,R 

= i[T L (e)f L (e) + r R (e)f R (e)]. (34) 
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Using the identities G r G a = {G r - G a )/{l/G a - l/G r ) = a(e)/T(e) [here o(e) = i[G r {e) 
G a (e)] is the spectral function], one can write G < in a 'pseudoequilibrium' form: 



G < {e) = ia{e)f{e) 



(35) 



where 



/» = 



T L (e)f L (e)+T R (e)f R (e 
Tie) 



a(e)= 



rfel 



eo-A(e)] 2 +[r(e)/2]2 " 
With these expressions the evaluation of the current (|T9|) is straightforward: 



J 



2^ 7 2tt 



a(6)[(r-(6)-r-(6))/(6) 

-(/L(e)r L (e)-^(e)r^)) 



(36) 



e r de 



r L (e)T R (e) 



h J 2 7r[e-e -A(e)] 2 +[r( e )/2] 5 



x[h(e) - f R (e)} . (37) 

Note that this derivation made no assumptions about proportionate coupling to the leads. 
The factor multiplying the difference of the fermi functions is the elastic transmission coeffi- 
cient. It is important to understand the difference between this result and the result obtained 
in Section III.B.2 (despite the similarity of appearance): There Eq . (|2l|) gives the current for 
a fully interacting system, and the evaluation of the retarded and advanced Green functions 
requires a consideration of interactions (e.g., electron-electron, electron-phonon, and spin- 
flip) in addition to tunneling back and forth to the contacts. Suppose now that the Green 
function for the interacting central region can be solved: G r,a (e) = [e — e — A(e) ±i r y(e)/2]~ 1 , 
where A and 7/2 are the real and imaginary parts of the self-energy (including interactions 
and tunneling). Then the interacting result for proportionate coupling ( pT|) becomes 
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e r de r „ , , . , x , r L (e)r*(e 



r^(e)+r«(e) 



[ e -6o-A( e )] 2 + [ 7 ( e )/2] 2 • 1 > 

This result coincides with the noninteracting current expression fl3"T|) if A(e) — > A(e) and 
7(e) — > T(e) = r /? (e) + T L (e). In a phenomenological model, where the total level width is 
expressed as a sum of elastic and inelastic widths, 7 = 7e + 7i, one recovers the results of 
Jonson and Grincwajg, and Weil and Vinter. pO 



C. Wide-band limit 

1. Basic formulae 

For simplicity, we continue to consider only a single level in the central region. As in the 
previous section, we assume that one can factorize the momentum and time-dependence of 
the tunneling coupling, but allow for different time- dependence for each barrier: Vk a (t) = 
UL/R(t)V atn (€k). Referring to Eq. fl33|) , the wide-band limit (WBL) consists of i) neglecting 
the level-shift A(e), ii) assuming that the line-widths are energy independent constants, 
J2a£L,R^a = r i//R , and (iii) allowing a single time-dependence, A L / R (t), for the energies in 
each lead. The retarded self-energy in Eq.(p9|) thus becomes 

57(M 2 )= £ <(ti)^(^)e" J ^ d ' 3Aa( ' 3) 

a£L,R 

x J ^-*to-*»)0(t 1 -t 2 )[-tT a ] 

= - i -[T L (t 1 ) + T R (t 1 )]5(ti-t2). (39) 

(Here we have introduced the notation Y L / R {t{) = r L / R (t 1 ,t 1 ) = r L / R \u L / R (ti)\ 2 .) With 
this self-energy, the retarded(advanced) Green function becomes |I7|J23]| 

G r > a (t, = g r > a (t, t') exp |=F jT dh 1 - [T^h) + T^h)}} (40) 
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with 

g r >*(t, t') = TiO(±t =f f) exp \-i f G?Wti)l . (41) 

L Jt' J 

This solution can now be used to evaluate the lesser function Eg. (|31~D , and further in Eg. (|15|). 
to obtain the time-dependent current. In the WBL the e- and ti-integrals in the term 
involving G < are readily evaluated, and we write the current as (using lm{G < (t, t)} — N(t), 
where N(t) is the occupation of the resonant level) 

JL(t)=- e j-{T L (t)N(t) + J^f L (e) 

x f dt 1 Y L {t 1 ,t)lnx{e^ t ^G r {t,t^}] . (42) 

J — oo ^ 

For a compact notation we introduce 

A L/R (e,t) = I dtiUL/nitjGrfah) 



x exp[te(t -h)-i f 1 dt 2 A L/R (t 2 )) . (43) 

J t 



N(t) = J2 rL/jR / 7^fL/n(e)\A L/R (e,t)\ 2 . (44) 



Obviously, in the time-independent case A(e) is just the Fourier transform of the retarded 
Green function G r (e). [|33| In terms of A(e,t) the occupation N(t) (using Eq.(|3~ID for G K ) is 
given by 

de 

L,R J 

We write the current as a sum of currents flowing out from the central region into the 
left (right) contact (see also Fig. ^|), and currents flowing into the central region from the 
left(right) contact, J L/R (t) = J%f R (t) + J% R (t): 

JTfk(t)= ~T L / R (t)N(t) (45) 
■#/«(*)= ~r L / R u L/R (t) J ^f L/R (e)lm{A L/R (e, t)} . (46) 



It is readily verified that these expressions coincide with Eq. (p7|) if all time-dependent guan- 
tities are replaced by constants. 
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Employing the same approach as in Section IV C, and provided that ttz,(t) = u R (t) = 

u(T), we find that the time-averaged current in the WBL is given by 

1e T L T R r de 

-f R (e)(u(t)A R (e,t))} . (47) 

Unlike the general case of Eq. fl27D , there is no restriction in the WBL that the energy 
dependence be the same in the two leads. Eq. ( [47|) can therefore be used for the case of 
a time-dependent bias, where A^(t) and A R (t) will be different. It is interesting to note 
that the function of energy appearing in the time-averaged current is positive definite. In 
particular, as is shown in Appendix C, 

- (lm{u L/R (t)A L/R (e,t)}) = l(\A L/R (e,t)\ 2 ) . (48) 

One consequence of (^Hj) is that if only the level is time dependent the average current cannot 
flow against the bias. 

In the next two sections we consider specific examples for the time variation, which are 
relevant for experimental situations. 



2. Response to harmonic modulation 

Harmonic time modulation is probably the most commonly encountered example of time 
dependence. Here we treat the case when the contact and site energy levels vary as 

A L /R,o(t) = A L/R>0 cos{ut) (49) 

It is easy to generalize the treatment to situations where the modulation frequencies and/or 
phases are different in different parts of the device. Assuming that the barrier heights do 
not depend on time {ul/r = 1), and substituting (|49| ) in the expression ( f4~3|) for A(e, t), one 
finds H 

. (A — Al/r) . 



A L /R(e, t) = exp [-i '- — sin(u;t)] 

CO 

ikuit 



10 

A - A L/R \ e 



k=— oo 



uj J e — eo — ku + iT /2 
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where J-k(x) = (— l) fc Jfc(x). Figures |3] (a) and (b) show \A(e, t)\ 2 and lmA(e,t) as a func- 
tion of time, respectively. We recall from Eqs.( fH| - f45|) that the current at a given time is 
determined by integrating \A(e, t)\ 2 and ImA(e, t) over energy, and thus an examination of 
Figure helps to understand to complicated time dependence discussed below. (We show 
only Al; similar results hold for Ar.) The physical parameters used to generate these plots 
are given in the figure caption. The three-dimensional plot (top part of figure) is projected 
down on a plane to yield a contour plot in order to help to visualize the time dependence. 
As expected, the time variation is periodic with period T = 2tt/uj. The time dependence is 
strikingly complex. The most easily recognized features are the maxima in the plot for \A\ 2 ; 
these are related to photon side-bands occuring at e = e ± ku (see also Eq.(|5T|) below). |35j 

The current is computed using the methods described in Appendix B, and is shown in 
Fig. f|. We also display the drive voltage as a broken line. Bearing in mind the complex 
time dependence of \A\ 2 and Iim4, which determine the out- and in-currents, respectively, 
it is not surprising that the current displays a non-adiabatic time dependence. The basic 
physical mechanism underlying the secondary maxima and minima in the current is the 
line-up of a photon-assisted resonant-tunneling peak with the contact chemical potentials. 
The rapid time variations are due to J m (or, equivalently, due to Iul4): the out-current J out 
is determined by the occupation N(t), and hence varies only on a time-scale T/h, which is 
the time scale for charge density changes. 

We next consider the time- average of the current. For the case of harmonic time depen- 



dence, we find [£34|] 



k=—oo 

X (e - e - ku) 2 + (r/2)2 ' (51) 
Figure [5] shows the resulting time-averaged current Jdc- A consequence of the complex 
harmonic structure of the time- dependent current is that for temperatures ksT < hu the 
average current oscillates as function of period 2tt/uj. The oscillation can be understood by 
examining the general expression for average current Eq . (|27D together with (|5lD: whenever 
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a photon-assisted peak in the effective density of states, occuring at e = eo ± koo in the time- 
averaged density of states (lmA L / R , moves in or out of the allowed energy range, determined 
by the difference of the contact occupation factors, a maximum (or minimum) in the average 
current results. 



3. Response to steplike modulation 

We give results for the case when the central site level changes abruptly at t — to: 
e o e o + A. If the contacts also change at the same time, the corresponding results are 
obtained by letting A — > A — A^m. Thus, simultaneous and equal shifts in the central 
region and the contacts have no effect. Assuming that the barrier heights do not depend on 
time {ul/r = 1), one finds for t > t from Eq.(f4~3"|) 

A(e,t) = ^— - 

e - e + u / 2 

l e - (e + A) + iV/2 J ' 1 ' 

This result is easily generalized (See Eq.(14) in Ref. [ |18|1 1) to a pulse of duration s, and 
numerical results are discussed below. 

It is instructive to study analytically the long and short-time behavior of A(e, t). It easily 
verified that A(e, t) has the expected limiting behavior: 

A(e,t -> oo) = [e- (e + A) +iY/2}- 1 . (53) 

Thus, when the transients have died away, A(e, t) settles to its new steady-state value. 

Consider next the change in current at short times after the pulse, t — 1 = 5t h/T, h/e. 
Note that the second inequality provides an effective cut-off for the energy integration re- 
quired for the current. In this limit we may write 

s 1 _ iA5t ,„ 

A(e,t)~ — — . (54) 

e - e + u 1 2 

Since 5J out (t) oc \A(e, t)\ 2 oc {St) 2 , the leading contribution comes from J m (t). For low 
temperatures we find 
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8J L/R (t) ~ — — / deIm6A(e,t) 
nn J-h/st 

~ — — Mt\og5t . (55) 
nn 

We next discuss the numerical results for a step-like modulation. Just like in the case 
of harmonic modulation, it is instructive to study the time dependence of \A\ 2 and Iul4; 
these are shown in Figures |](a) and (b), respectively. The observed time dependence is 
less complex than in the harmonic case. Nevertheless, the resulting current, which we have 
computed for a pulse of duration s, and display in Fig. |7], shows an interesting ringing 
behavior. The ringing is again due to the movement of the sidebands of 1mA l/r through 
the contact Fermi energies. 

Due to the experimental caveats discussed in Section II, the ringing showed in Fig. [7] 
may be masked by capacitive effects not included in the present work. However, the ringing 
should be observable in the time-averaged current by applying a series of pulses such as that 
of Fig. [5], and then varying the pulse duration. |36|] In Fig. [8] the derivative of the dc current 
with respect to pulse length is plotted, normalized by the repeat time r between pulses. For 
pulse lengths s of the order of the resonance liftime %/Y, the derivative of the dc current 
mimics closely the time-dependent current following the pulse, and, likewise, asymptotes to 
the steady-state current at the new voltage. 

4- Linear-response 

For circuit modeling purposes it would often be desirable to replace the mesoscopic 
device with a conventional circuit element, with an associated complex impedance Z(u), or 
admittance Y(uS). Our results for the nonlinear time-dependent current form a very practical 
starting point for such a calculation. For the non-interacting case, the current is determined 
by A(e,t) (see Eq.( ^ - |45|) ), and all one has to do is to linearize A (Eq. fl43[) ) with respect to 
the amplitude of the drive signal, i.e., A — Al/r- It is important to notice that we do not 
linearize with respect to the chemical potential difference: the results given below apply to 



25 



an arbitrary static bias voltage. 

Performing the linearization, one finds 



\A 



(i) 



e,t)| 2 = ^^Re{ ^ 

' ' (,) I- f — -+- 7.T 



]}, (56) 



and 



e — eo — uj — iT/2 e — to + uj — iT/2 



lmA {1) (e t) - A ~ AL/i W e "~" 
LmA L/R (e,t)- ^ M c _ ^ _ w + ir /2 

~e-e + u + iT/2 + e-e + iT/2^ ' ^ 



At finite temperature the energy integration must be done numerically, as explained in 
Appendix B, while at T = they can be done analytically. In the latter case, all the 
integrals can be cast into the form 

de 



f — 

J-oo (e — 



(e-ei + iT 1 /2){e-e 2 + iT 2 /2) 



log "- ei+i ^ . (58) 



ei-e 2 + i(r 2 -ri)/2 & /i - e 2 + iY 2 /2 
Using log(x + iy) = l/21og(:r 2 + y 2 ) + itan^^y/x) yields 



+ sin(wi)G ! L/ fl(w)] (59) 



and 



j(l),out _ e pL/fi pL/il ^ ^L/J/ 



{cOS^I ^^^ Cl/rH - ^^p ^L/flM] 

r uj 

— sm(ut) [ — — ^ 2 Gl/r{uj) + ^ 2 ^ r2 F £/jR (o;)]}, 



(60) 



where we defined 
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Gl/r{u) = log ' ^ 61 



and 

-i ^l/r - e Q — u 



Fl/r,{u) = tan" 



T/2 



-tan- 1 ^-^o + o; 

r/2 v ; 

These expressions give the linear ac current for an arbitrarily biased double barrier structure, 
where both contacts and the central region energies are allowed to vary harmonically. As a 
check, it is instructive to verify that the finite temperature results of Appendix B.2 contain 
Eqs.( [59| - |60| ) as a special case; this is a rather straightforward calculation using the limiting 
behavior of the Digamma function. 

Considerable simplification occurs, if one considers a symmetric structure at zero bias: 
r L = T R = T/2, and fiL = fiR = A*, respectively. Following Fig. |], the net current from left 
to right is J« = 1/2 [4 1)lin + 4 1} ' out - 4 1)>out - 4 1),in ]. Using Eq.(]5j|fD, one finds that the 



'out' contributions cancel, and that the 'in'-currents combine to give the net current 

J«(t) = -§£ Al ~ AR [cos(ojt)F(uj) + wy(ut)G(u)] . (63) 
n 4 2ttuj 



Here the functions F(uj) and G(u>) are given by Eq. ( |5"2"| - |51~D , but using \i and r/2 as param- 



eters. This result exactly coincides with the recent calculation of Fu and Dudley, p2| which 
employed a different method. 

We now wish to apply the formal results derived in this section to an experimentally 
relevant system. The archetypal mesoscopic with potential for applications is the resonant- 
tunneling diode. The key feature of a resonant-tunneling diode is its ability to show negative 
differential resistance (NDR). The WBL model studied in this section does not have this 
feature: its IV-characteristic, which is readily evaluated with Eq. (|37D , is a monotonically 
increasing function. A much more interesting model can be constructed by considering a 
model where the contacts have a finite occupied bandwidth; this can achieved by introducing 
a low energy cut-off Dl/r (in addition to the upper cut-off provided by the electro-chemical 
potential). The zero-temperature IV-characteristic is now 
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Acy ' h r L r/2 r/2 

/Xfl(V) - e (VQ _ t ti R (V)-D R -e (V) , 

-tan ^ +tan ^ j. 

(64) 

Here we assume that the right chemical potential is field dependent: Hr{V) = Hr — eV, 
and that the field-dependence of the central region level is given by e (V) = e — V/2. 



The resulting current-voltage characteristic is depicted in Fig. |10[ We note that the strong 
increase in current, which is observed in experimental systems at very high voltages, is not 
present in our model: this is because we have ignored the bias-dependence of the barrier 
heights as well as any higher lying resonances. The only generalization required for Eqns. (|59| 



-50|) is to modify the F and G functions: F^ —> F = F^ — F^ D , and analogously for G^. 
We show in Fig. |ll] the resulting linear-response admittance Y(u) for a symmetric structure 
(T^ = Tr). Several points are worth noticing. For dc bias eV = 5 (energies are given in 
units of T) the calculated admittance resembles qualitatively the results reported by Fu and 
Dudley for zero external bias, except that the change in sign for the imaginary part of Y(u) 
is not seen. For zero external bias (not shown in the figure) our finite band-width model 
leads to an admittance, whose imaginary part changes sign, and thus the behavior found by 
Fu and Dudley cannot be ascribed to an artefact of their infinite band-width model. More 
interestingly, for dc bias in the NDR regime, the real part is negative for small frequencies. 
This simply reflects the fact that the device is operating under NDR bias conditions. At 
higher frequencies the real part becomes positive, thus indicating that further modeling along 
the lines sketched here may lead to important implications on the high-frequency response 
of resonant-tunneling structures. 

In concluding this section, we wish to emphasize that the linear-response analysis pre- 
sented above is only a special case of the general results of Section IV, which seem to have 
the potential for many applications. 
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VI. RESONANT TUNNELING WITH ELECTRON-PHONON INTERACTIONS 



As a final application, we establish a connection to previous calculations on the effect of 



phonons on resonant tunneling. [23] For simplicity, we consider a single resonant level with 
energy-independent level widths Tl and Tr (i.e. the WBL). The expression for the current 
Eq.(PTD becomes now 

where a(t) = i[G r (t) — G a (t)} is the interacting spectral density. In general, an exact eval- 
uation of a(t) is not possible, however, if one ignores the Fermi sea, G r (t) (and hence a(t)) 
can be calculated exactly: f3q| 

G r (t) = -i9{t) exp[-it(e - A) - - Tt/2] , (66) 



where 



M 2 

A = E— > (67) 



q ^q 



and 



M 2 

$ W = E "^[^q(l - e " q< ) + (^q + !)(!- e ^ q *)] > (68) 



,2 

q ^q 



and the electron-phonon interaction is given by Eq.(||). When substituted in the expression 
for current, one recovers the result of Ref. [ |23|]] , which originally was derived by analyzing the 
much more complex two-particle Green function G(t, s, t) = 9(s)9(t)(d{r— s)^ (rfdfyS (0)) . 
The advantage of the method presented here is that one only needs the smg/e-particle Green 
function to use the interacting current formula (|2l|). Other systematic approaches to the 
single-particle Green function can therefore be directly applied to the current (e.g. pertur- 
bation theory in the tunneling Hamiltonian) . 
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VII. CONCLUSIONS 



Here, we summarize the main results of this study. We have derived a general formula for 
the time-dependent current through an interacting mesoscopic region, Eq. flT5|) . The current 
is written in terms of local Green functions. This general expression is then examined in 
several special cases: (i) It is shown how earlier results for time-independent current are 
contained in it [Eqs. (|I^) and (Ell)], (ii) An exact solution, for arbitrary time- dependence, 
for a single non- interacting level coupled to two leads is given [Eqs. (E3j - E5|) ]. This calculation 
leads to a prediction of 'ringing' of current in response to abrupt change of bias, or in response 
to an ac-bias. We believe that this prediction should be experimentally verifiable, (iii) We 
derive a Landauer-like formula for the average current, Eq. fl27|) , and apply it to a simple 
model of dynamical disorder. Finally, as applications, we discuss (iv) a.c. linear-response 
at arbitrary dc-bias and finite temperature, and (v) find a connection to earlier results on 
resonant tunneling in the presence of optical phonons. 

We hope that time-dependence will provide a new window on coherent quantum- 
transport, and will lead to significant new insights in the future. 
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APPENDIX A: NONEQUILIBRIUM GREEN FUNCTIONS 



The most important result (see, e.g., Refs.[ |2"5|] , [CT, P7[]) of the formal theory of 
nonequilibrium Green functions is that the perturbation expansion has precisely the same 
structure as the T = equilibrium expansion. Instead of a time-ordered Green function, 
one works with the contour-ordered Green function, 

G(T,T') = -i(T C {lP(T)^(T')}), (Al) 

where the contour C is shown in Fig. 0. The contour-ordering operator Tc orders the 
operators following it in the contour sense: operators with time labels later on the contour 
are moved left of operators of earlier time labels. Thus, once the self-energy functional, 
S = £[£?], has been specified, the contour-ordered Green function obeys formally the same 
Dyson equation as in T = theory, 

G = Go + G SG , (A2) 

with the modification that internal time-integrations run along the (complex) path discussed 
in section II. A. It follows from this structural equivalence that one can derive equations of 
motion just as in the T = case, and that the passage to nonequilibrium takes place by 
replacing the time-ordered Green functions by contour-ordered Green functions, and by 
replacing the real-time integration by an integration along the time-contour. In practical 
calculations, however, the contour ordered Green functions are inconvenient, and it is expe- 
dient to perform an analytic continuation to the real axis. The first step in this procedure 
consists of expressing the contour ordered Green functions in terms of 2 x 2 matrices, whose 
elements are determined by which branches of the contour the two time labels are located 
on. The four elements of the matrix Green function are not linearly independent, and it is 
useful to perform a rotation of this matrix. A particularly convenient set of operational rules 



has been given by Langreth: [2^] If one has an expression A — J BC on the contour (this is 
the generic type of term encountered in the perturbation expansion), then the retarded and 
lesser components are given by 
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A r (t,t')= J dt 1 B r ^,t x )C r ^x,H) 
A<(t,t')= J dt^^t^^f) 

+B<{t,t 1 )C a {t 1 ,t')) (A3) 

These results are readily generalized to products involving three (or more) Green functions 
or self-energies. 

The equation of motion for G K can be derived by applying the rules ( |A3| ) to the Dyson 
equation for the contour-ordered Green function. The Dyson equation can be written either 
in a differential form, or in an integral form, as in Eq . flA"2]) . The former leads to the Baym- 
Kadanoff transport equation, while the latter (which is employed in the present work) yields 
the Keldysh equation for the lesser function: 

G< = (1 + G r S r )G<(l + E a G a ) + G r S<G a , (A4) 

where the retarded and advanced Green functions satisfy 

G r,a = G r,a + Q^r^a _ ^ 

The physical modeling goes in the choice of the self-energy functional E, which contains the 
interactions (carrier-impurity scattering, phonon scattering, carrier-carrier scattering etc.). 
Once E is given, for example in terms of diagrams, the retarded, or 'lesser' components of 
the self-energy can be worked out according to the rules (|A3|) , and one can proceed to solve 
the coupled equations (|A4|) and (|A5|) . 



APPENDIX B: DYSON EQUATION FOR G t NJ(a 

1. Equation-of-motion method 

According to Appendix A it is sufficient to consider the T = equation of motion for 
the time-ordered Green function G^ fcQ : 
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9 it 

~ ^g~^G nka (t ~ t ) = e kG n ka (t - t ) 

m 

where we defined the central region time-ordered Green function function G^ m (£ — t') = 
— i(T{dJ n (t')d„(t)}). Note that it is crucial that the leads be non-interacting: had we 
allowed interactions in the leads the equation of motion technique would have generated 
higher order Green functions in Eq . ( [BT|) , and we would not have a closed set of equations. 

We can interpret the factors multiplying G\ ka (t — t') as the inverse of the contact Green 
function operator, and introduce a short-hand notation: G^^g^ = J2 m G t nm V kam . By 
operating with g\ a from right, we arrive at 

G^kait - 1') — 22, I dtiGr nm (t - h ) V£ a >m g ha (t 1 - t') . (B2) 

m 

According to the rules of the nonequilibrium theory, this equation has in nonequilibrium 
precisely the same form, except that the intermediate time integration runs on the complex 
contour: 

G nM (T,T') =J2j rfr i G nm(r,ri)\4* aim (ri)^ fea (ri,r') . (B3) 



This is Eq . (|TT|) of the main text. The analytic continuation rules (|A3|) can be applied, and 
the desired Dyson equation is obtained. 

2. 5-matrix expansion 

We write the Green function G n ^ a {t, t') in terms of interaction-picture operators (denoted 
by a tilde) by invoking the S'-matrix: 

G nM (r, r') = -t(T c {Sd n (r)~c{ a (r')}} , (B4) 

where 

S = T c {exp[-i [ dnHrin)}} (B5) 

J C 
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is the contour-ordered S'-matrix, and Ht is the tunneling Hamiltonian of Section II. B2. We 
expand the exponential function in (|B~5|); the zeroth order term does not contribute, and we 
find 



G nM (T,T') = -i{T c [d n {T)c\ a {r') J2 



00 < —j\n+l 



X 



[/ dT< 2 \ V kiai, m {T2)c\, a ,{T 2 )d m {T2) 
^ k'a',m 

+nV,m(T2)4(r2)c fc w(r 2 )]] n+1 }) . (B6) 



Since, by assumption, the leads are non-interacting, result will only be non-zero if c\. a (r') is 
contracted with Cfc Q (rj) from one of the n + 1 interaction terms. The n + 1 possible choices 
cancels a factor of n + 1 in the factorial in the denominator, leaving 

G nM (r,r')=J2 f dr 2 (-z)(T c {c fcQ (r 2 )cL(r')}) 

xVL m (r 2 )(-i)(T c {Sdl(T 2 )l(r)}) . (B7) 



Eq.( p7|) is completely equivalent to the result (|B3|) obtained in the previous subsection. 



APPENDIX C: PROOF OF EQ.(pD 



In this Appendix, we prove that for a single level in the WBL (see Section V C) there is 
a definite relation, 

- (u L/R (t)lm{A L/R (e,t)}) = ^(\A L/R (e,t)\ 2 ), (CI) 

between the time averages of the quantities that, respectively, determine the current and 
the occupation. For the case of the occupation, one can explicitly write out (\A L / R (e, t)\ 2 ) 
and then use the identity 

G r {t, ti)G°(4 1) = i6{t - t x )9{t - t[) 

x 'e-W-^Grfa, h) - e- r( *-* l) G a (t' 1 , h)] (C2) 

to obtain 
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2 rl /Z rl /Z 

(|A| 2 ) = hm — y_ T/2 ^i ]_ T/2 d1& L , R {h)u L / R {&) 
x{G r {t' x M) ~ G a {t'iM)) exp[ie(«i - ti) + dt 2 A(t 2 )]. 



(C3) 



Writing out (w£,/#(£)Im{A£/#(e, i)}) explicitly then yields Eq. (|(J1[). 



APPENDIX D: NUMERICAL INTEGRATION 

In this Appendix we describe methods to facilitate numerical calculations in the wide- 
band limit (Section |VC ). While the numerical integrations required for the occupation and 



for the current can be done directly, it is often difficult to obtain sufficient accuracy. We 
have found that it is useful to do the integrations analytically by contour integration, and 
then sum the resulting residues. We have also checked for a few selected parameter values 
that the two methods give identical results. 

1. Steplike modulation 

We illustrate the somewhat cumbersome but straightforward formulae by giving the 
expressions for the deviation of the occupation from its asymptotic value following a steplike 
modulation of the level energy (Section |VC3|) : 5N(t) = N(t) — N(t = oo). We find from 
Eqs.flU) and (0) 

SN(t) = iAV r H)[r^(^) + r fl Dy] 

271 

_±Ae- r ^)/ 2 [r L 2Re{E(/i L )} + T R 2Bje{E(fi R )}] , (Df ) 

Z7T 



where 



D(fi) = J de- { 
E(n) = J de[- 



f(e) 



e - e - A)^ + (172)2 (e - e o y + (r/2) 



2 



J{e-e -A)(t-t ) 



e - e - A)2 + (r/2)2 e - e - ^/2 
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(D2) 



where /(e) is the Fermi function with chemical potential \x. The poles of the integrands are at 
e = e ±iT/2, e = e + A±ir/2, and e = fi±i2n(n + 1/2)//?, respectively. Upon closing the 
contour in the upper-half-plane, one obtains three different contributions; the terms arising 
from e = €q + iT /2 and e = eo + A + z7/2 obviously lead to no problems, while the sum over 
n converges either as n~ A (the term originating from D(fi)), or as n~ 3 exp(— 2im(t — to)/f3) 
(the term due to E(/j)), and hence also converges rapidly. 



2. Harmonic modulation 



In principle, the calculation proceeds as in the previous section. However, the sum over 
the residues, which results from the contour integration, converges very slowly. A typical 
term in the resulting lengthy expressions converges only as n~ 2 . Significantly improved 
convergence can be obtained by making use of the relation 

OO 1 1 

y 1 = _ 

^ {n + a)(n + b) a- 



■[*(a) - * 



(D3) 



where \& is the digamma function. In what follows, we give the results for linear-response. 



The occupation [which also gives the current flowing out from the central region via fl44j) 1 is 



1 „A -A z/fl r £ /* r 



L.R 

r L/R _ jL/R 



10 



sin(cut)(2rro /R 



1 

+u>{lZ' K - It'") - T(R L J R + R L _ /R )) 
+ cos{ujt)(-2ujr R/L + uo{R L J R + R R,L ) 



+T{1 



L/R jL/R 



)) 



(D4) 



Here 



I L ± /R = Im 
Rl /R = Re 



Tq — Re 



(1/2 - -^-Xhl/r - e T w ~ iT/2)) 
^(1/2 - -^zip-L/R -to Too- ir/2)) 
*(l/2 + ^(fiL/R-eo + iT/2)) 



(D5) 



The current flowing into the central region can also be expressed in terms of similar functions: 
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with 



■#/*(*)= ^ L/R 



.A - A 



L/R 



llXOJ 



cos(ujt)(i L ^ R — i+ R ) 



+ sin(ujt)(2rQ^ R — r+ R — r^ R )\ , 



(D6) 



i L J R = Im[*(l/2 + ^(^L/i? - e =F + ir/2)) 
r^ /i? = Re[*(l/2 + /r(/^/* ~e TuJ + ir/2))' 



(D7) 



By recalling lim^^ ^(z) — > log(z), it is straightforward to check that these results reduce 
to the T = case discussed in the main text. 
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FIGURES 

FIG. 1. Sketch of charge distribution in a three-dimensional resonant-tunneling device under 
dc-bias Vbias = PL ~~ PR with a time-modulation of amplitude £±l/r superposed on the leads. As 
argued in the text, only a tiny fraction of charge carriers participates in setting up the voltage drop 
across the structure. 

FIG. 2. The complex-time contour on which nonequilibrium Green function theory is con- 
structed. In the contour sense, the time t\ is earlier than T2 even though its real time projection 
appears larger. 

FIG. 3. (a)|j4(e, t)\ 2 as a function of time for harmonic modulation for a symmetric structure, 
Tl = = T/2. The unit for the time-axis is ti/T, and all energies are measured in units of T, 
with the values hl = 10, hr = 0, eo = 5, A = 5, A^ = 10, and A# = 0. The modulation frequency 
is oj = 2F/H. (b) The time-dependence of lmA(e,i) for the case shown in Fig 3(a). 

FIG. 4. The time-dependent current J(t) for harmonic modulation corresponding to the pa- 
rameters of Figure 3. The dc bias is defined via hl = 10 and = 0, respectively. The dotted 
line shows (not drawn to scale) the time dependence of the drive signal. The temperature is 

k B T = o.ir. 

FIG. 5. Time averaged current J^ c as function of the ac oscillation period 2tt/uj. The dc 
amplitudes are the same as those in Fig. 4. 

FIG. 6. (a) \A(e, t)\ 2 as a function of time for step-like modulation. At t = the resonant-level 
energy eo suddenly decreases by 5r. (b) The time dependence of Iim4(e, t) for the case shown in 
Fig 6(a). 
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FIG. 7. Time-dependent current J(t) through a symmetric double-barrier tunneling structure 
in response to a rectangular bias pulse. Initially, the chemical potentials hl and [1r and the 
resonant-level energy eo are all zero. At t = 0, a bias pulse (dashed curve) suddenly increases 
energies in the left lead by Al = 10 and increases the resonant-level energy by A = 5. At t = 3, 
before the current has settled to a new steady value, the pulse ends and the current decays back 
to zero. The temperature is ksT = 0.1T. 

FIG. 8. Derivative of the integrated dc current J,ic with respect to pulse duration s, normalized 
by the interval between pulses r. For pulse durations much longer than the resonance lifetime h/T, 
the derivative is just the steady-state current at the bias voltage, but for shorter pulses the ringing 
response of the current is evident. 

FIG. 9. Linear-response configuration. 

FIG. 10. IV-characteristic for a model resonant-tunneling device (quantum dot). The system 
is defined by parameters eo(V = 0) = 2, //£ = hr{V = 0) = 0, and Dl = Dr = 2, and the current 
is given in units of eT/h. 

FIG. 11. In-phase and out-of-phase components of the linear response current (in units of eT/h 
and normalized with the amplitude of the drive signal to yield admittance) for two bias points, 
eV = 5 (continuous line) and eV = 10 (dashed line). Other parameters are as in Figure 10. 
The out-of-phase components (or, equivalently, imaginary parts) always tend to zero for vanishing 
frequency, while the in-phase component can have either a positive or negative zero-frequency limit 
depending on the dc bias. 
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